Time-reversal symmetry breaking phase in the Hubbard model: a VCA study 
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We study the stability of the time-reversal symmetry breaking staggered flux phase of a single band 
Hubbard model, within the Variational Cluster Approach (VCA). For intermediate and small values 
of the interaction U, we find metastable solutions for the staggered flux phase, with a maximum 
current per bond at U ~ 3.2. However, allowing for antiferromagnetic and superconducting long- 
range order it turns out that in the region at and close to half filling the antiferromagnetic phase 
is the most favorable energetically. The effect of nearest-neighbour interaction is also considered. 
Our results show that a negative nearest-neighbour interaction and finite doping favors the stability 
of the staggered-flux phase. We also present preliminary results for the three-band Hubbard model 
obtained with a restricted set of variational parameters. For this case, no spontaneous time-reversal 
symmetry breaking phase is found in our calculations. 
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I. INTRODUCTION 



The Hubbard Hamiltonian is perhaps the simplest pos- 
sible model to describe electronic interactions in solids. 
It is however, difficult to solve and even after several 
decades of research there are still many open questions 
about its basic features. In particular many competing 
phases have been proposed as candidates for the ground 
state in special parameter regions. The staggered flux 
phase [H, 0], known also as "orbital antiferromagnet" Q, 
was proposed to describe the ground state of Heisen- 
berg Antiferromagnetism. Many properties of the stag- 
gered flux phase have been discussed in the past. How- 
ever, less numerical evidence exists that support the ex- 
istence of this phase in the Hubbard or t-J model 0- 
[To| . The phenomena of pseudogap in cuprate supercon- 
ductor contributed to the increased interest in the time- 
reversal symmetry broken phase. Chakravarty et al. 
argued that the pseudogap region is due to the compe- 
tition between d-wave superconductivity and d-density 
wave (DDW) state [HBj, which actually is the stag- 
gered flux state breaking both time-reversal and trans- 
lational symmetry. Based on the mean-field analysis of 
three-band Hubbard model, Varma showed that a time- 
reversal symmetry breaking phase, the "circulating cur- 
rent phase" [14H17I ]. is stable in some parameter regions. 
The orbital current of this phase is circulating along the 
O-Cu-0 plaqette in each cell and thus, in contrast to 
DDW, translational symmetry is preserved. Motivated 
by this proposal, several theoretical and experimental 
groups have tried to find the signatures of orbital cur- 
rent in the CuO planes [l8l - l2l| . However, no agreements 
have been obtained so far. Recent polarized neutron 
diffraction experiment shows that the spontaneous cur- 
rent occurs in loops involving apical oxygen orbitals [22j . 
This picture is consistent with Variational Monte Carlo 



(VMC) calculations [23[, although the authors observe 
that the computed current of three-band Hubbard model 
is smaller, for larger system sizes, therefore computations 
on larger system are required. 

There have been a lot of efforts to search for the stag- 
gered flux phase in the two-leg ladder, which is easier 
to study and can shed some light on the full two dimen- 
sion lattices. By using the highly accurate density-matrix 
renormalization group technique, Marston et. al. [24| 
and Schollwock et al. [25[ found the evidences for the 
existence of staggered flux phase at and away from half 
filling, respectively. It has been shown that complex in- 
teractions are needed to stabilize this phase [24|,[25j]. Ana- 
lytical studies using bosonization/renormalization group 
(RG) method also found stable regions of this phase 
for weak interactions 26- 29]. In the present work we 
adopt the Variational Cluster Approach (VCA) to study 
the time-reversal symmetry breaking phase. VCA allows 
to evaluate the single-particle [3(| and two-particle [3lJ 
spectral functions, and, to include symmetry breaking 
phases [3(| HH, HH , such as Antiferromagnetism and Su- 
perconductivity, therefore it can be easily extended to 
include time-reversal symmetry breaking phase. 

The paper is organized as follows. In Sec. [Hi we 
present the Hubbard Hamiltonian in a form that includes 
the single particle parameters that are relevant for the 
time-reversal symmetry breaking. We review briefly the 
VCA method and introduce the coupling fields associ- 
ated with the time-reversal symmetry breaking, so called 
"Weiss" fields. These are the hopping At and its phase 
A(j>, together with the current per bond which arises as 
a natural order parameter. The search for the sponta- 
neous time-reversal breaking state within the interact- 
ing Hubbard model is presented in Sec. IIIII The follow- 
ing sections Sees. IIII Al and MI Bl discuss the stability 
of staggered flux phase at and away from half-filling for 
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the single band Hubbard model. Our results indicate 
that the stable ground state for the interacting Hubbard 
model at half filling is the antiferromagnetic state. Away 
from half-filling, the competition with the superconduct- 
ing phase is also investigated. Sec. IIII CI discusses the 
effect of nearest-neighbour interaction on the staggered- 
fiux phase. In Sec. IIII D[ we extend the procedure de- 
scribed in Sec. [IT] to the three band Hubbard model and 
search for a spontaneous time-reversal symmetry break- 
ing phase in this model. Finally, we draw our conclusions 
in Sec. UVl 



II. HAMILTONIAN AND VARIATIONAL 
CLUSTER APPROACH 

We consider the following 2D Hubbard model on a 
square lattice: 



H 



<ij>cr 



+ h.c.) + U riifnn (1) 



where cj CT (ci CT ) creates (destroys) an electron on site 



with spin er, 



is the particle number opera- 



tor, (ij) denotes nearest neighbor sites, and U is the local 
Coulomb interaction, tij — te l ^ ij is the hopping-matrix 
element, where t is a real number and <pij is the phase 
factor on the bond (ij). In the following, we take t = 1 
to be our unit of energy. If all fyj = 0, Eq. (UJ is the 
Hamiltonian of the standard Hubbard model with real 
hopping matrix. For <^>y ^ 0, Eq. (TTJ) describes electrons 
hopping on lattices subjected to a magnetic field perpen- 
dicular to the plane of lattices provided the circulation of 
4>ij is nonzero. For <pij = the Hamiltonian is symmetric 
under time reversal. 

The Variational Cluster Approach (VCA) [H HI em- 
ployed in this paper can be seen as an extension of Cluster 
Perturbation Theory (CPT) [Hr{37| based on the Self- 
energy Functional Approach (SFA) [38|, [3{|. One short- 
coming of CPT method is that it can not study the spon- 
taneous symmetry-breaking phase. As an improvement 
over it, the VCA adds additional "Weiss" field (describ- 
ing a particular ordered state) to the cluster Hamiltonian, 
and optimizes the self-energy of the reference cluster by 
variational principle. The value for the variational pa- 
rameter is determined by searching the saddle point of 
SFA grand potential, 

fi = Q' + Trh^Go" 1 - S) _1 - TrlnG' (2) 

where Go is the non- interacting Green's function, 57', X, 
and G' are the grand-canonical potential, self energy, and 
Green's function of the reference cluster, respectively. 
The system is in the ordered state if the corresponding 
"Weiss" field is nonzero at the saddle point of grand po- 
tential n. 

For an initially time-reversal symmetric Hamiltonian 
(4>ij — 0) , spontaneous symmetry breaking can occur un- 
der some circumstances. Within VCA this can be stud- 
ied by introducing in the reference system an appropriate 



"Weiss" field [30|]. The "Weiss" field (H^ R ) that breaks 
the time-reversal symmetry of the cluster Hamiltonian 
can be written in the following form: 



E 

<ij>a 



(At e tA ^clc Jr7 + h.c.) 



(3) 



Here, At is the amplitude of hopping acting as the 
strength of the "Weiss" field, and Afcj is the phase factor 
on bond (ij) which determines the phase configuration 
of time-reversal symmetry-breaking phase. Both At and 
A(f>ij should be treated as variational parameters when 
searching the saddle points of the SFA grand potential 
in Eq. @. A time-reversal symmetry-breaking phase is 
characterized by a saddle point of f2 with nonzero values 
of both At and A(j>ij . 

The natural order parameter for the time-reversal 
symmetry- breaking phase is the current. For a typical 
bond (lm), the current is defined as 



Jim = ^ it(< 4 



> - < ci 



7 Cla >) 



(4) 



It is obvious that Ji m is nonzero only if the time-reversal 
symmetry is broken. In the staggered flux state the phase 
displays a pattern like the one illustrated in Fig. [TJ If no 
other symmetries are broken one expects the absolute 
value of the phase to be the same on all bonds. We thus 
take |A0jj| = A<fi to be the same on all bonds [40j . 

Before proceeding to search for the spontaneous time- 
reversal symmetry breaking phase, we first test our code 
by introducing an external staggered flux phase on a 
square lattice induced by an external magnetic field. The 
dynamics of two dimensional lattice electrons coupled 
to an external magnetic field is described by the tight- 
binding Hamiltonian in Eq. (UJ), with the circulation of 
4>ij on a given path being proportional to the magnetic 
flux through the enclosed surface [4l| . 
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FIG. 1: (Color online). Staggered flux pattern in a 2D square 
lattice. <f> is the phase on each bond in the direction indicated 
by the arrow. The total magnetic flux through a plaquette 
is 40. + and — denote the directions of the corresponding 
magnetic field. The black dots denote the 10-site reference 
cluster we use in our calculations. The vectors e*i and e*2 are 
the translation vectors of VCA. 
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For a staggered flux pattern and interaction [7 = 0, the 
model can be solved analytically and the energy spec- 
trum obtained. The ground state Eq and current J 
can be computed by summing all the energies up to the 
Fermi energy and by the Hellmann-Feynman theorem. 
We carry out our test in the non-interacting limit for a 
staggered magnetic field with a phase pattern as shown 
in Fig. [TJ All the quantities evaluated: spectral function, 
Eg{4>), and J using our VCA code agree precisely with 
the analytical results. 

For the case of interacting electrons we additionally 
checked the validity of the implementation by verifying 
thermodynamic consistency. Specifically concerning the 
current, its expectation values evaluated via the spec- 
tral function and via the Hellmann-Feynmann theorem 
turn out to be identical. Although this is obviously ex- 
pected, it in practice non trivial since, as discussed in 
42] , this is only guaranteed if one uses the correspond- 
ing hoppings as variational parameters. By separating 
the hopping term in Hamiltonian (fTJ into real and imag- 
inary parts: ct(c\ a c ]a + ct a c il7 ) and i^ ia c ja - J ja c ia ), 
with a = — icos(0) and j3 — — tsm(<^>), we confirm that 
treating j3 as variational parameter the thermodynamic 
consistency condition discussed above is satisfied. 



III. SPONTANEOUS STAGGERED FLUX 
PHASE IN THE 2D HUBBARD MODEL 




FIG. 2: (Color online). Difference between the grand poten- 
tials computed with variational parameters, such as the hop- 
ping At, phase factor Acj>, and staggered magnetic field Iiaf, 
with respect to the grand potential Q(0) when no variational 
parameters are considered. The calculations are performed at 
half filling. 

In this section, we study the stability of the sponta- 
neous staggered flux phase in the Hubbard model with 
real hopping matrix (i.e., <f>ij = in Eq. ([I])). The situa- 
tion is different from what we discuss in Section [TTI where 
the electrons are under an external staggered magnetic 



field. In order to allow for the time-reversal symmetry- 
breaking order, we introduce a Weiss field H' TR (cf. Eq. 
©) in the reference system. For the staggered flux phase, 
the absolute value of the phase |A<^y| is chosen to be 
equal for all bonds (i.e, A</>) and the staggering configu- 
ration is shown together with the VCA reference cluster 
in Fig. [TJ We stress that it is important to choose a 
reference cluster with even number of plaquettes, so that 
the net "magnetic field" in the cluster is zero. 

A. Staggered flux phase in the Hubbard model at 
half-filling 

In Fig. [5] we plot the difference between the values 
of the grand potential obtained when variational param- 
eters are included relative to the cases when no varia- 
tional parameters are considered. The results presented 
are for the case of half-filling and different values of U. 
In the small U region (i.e., from to 3.6), the values 
of the grand potential S7(At, A<p) with both At and A<j> 
taken as variational parameters are always smaller than 
the grand potential SI (At) with a single variational pa- 
rameter At, the latter being also smaller than the grand 
potential obtained in a calculation with no variational 
parameters S1(0). Accordingly, the staggered flux phase 
is lower in energy than the correlated paramagnetic (non- 
magnetic) state of the tight-binding Hamiltonian at half- 
filling. At U ~ 3.6 the optimal A(j> vanishes indicating 
a phase transition to a phase without time-reversal sym- 
metry breaking. This suggests that the staggered flux 
phase is favored by small values of U. This situation was 
not considered within Ref. [1], which predicts that the 
staggered flux phase is stable at large U limit, a conse- 
quence of the mean-field approach. At half filling one 
should expect that the ordered Antiferromagnetic phase 
also plays an important role in the Hubbard model. In or- 
der to compare the stability of the flux phase with respect 
to the Antiferromagnetic long-range order we included a 
Neel Antiferromagnetic Weiss field 

H'af = h AF £> <t - n 4 )e^ R < (5) 

i 

into the reference system. Here, Iiaf is the strength of 
a staggered magnetic field and Q = (tt,tt) is the Anti- 
ferromagnetic wave vector. The calculations within the 
long-range ordered state shows that SI (At, A<p) is larger 
than fl(hAF) for values of U > 1.9, and is lower than 
0(/!af) when U < 1.9. But for a complete comparison, 
fl(At, Acf>) has to be compared with f2(Ai, Haf) when 
both At and Kaf are treated as variational parameters 
in the long-range ordered state. As shown in Fig. [21 in 
this case SI (At, A<\>) is always larger than Q(At, Haf) for 
the whole range of U. 

Although our results indicate that the Antiferromag- 
netic state is more stable than the staggered flux phase 
for the half-filled 2D Hubbard model, the latter can be 
still considered as a metastable phase whose fluctuations 
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affect the physics of the system, or which may become 
stable whenever Antifcrromagnetic long-range order is 
suppressed by some other mechanism. 

In order to study the metastable staggered flux state 
in more detail, we plot in Fig. [3ji the difference in grand 
potentials i}(At,A<f>) — fl(At) (symbol - □, with label- 
ing on the left) and the corresponding current J as a 
function of interaction U at half filling. The difference 
between the values of the grand potential increases with 
U, reaches a maximum for the value of U « 3.2, and then 
decreases to zero for U ~ 3.6. The largest difference be- 
tween fi(At, A(j>) and fl(At) is 0.0027 which corresponds 
to about 1 meV for t w 500meV. The current per bond 
J (obtained in the staggered flux phase) is computed by 
treating both At and A(f> as variational parameters. The 
computed values are shown in Fig. [3^ with labeling on 
the right. One sees an opposite behaviour of the current 
as a function of U with respect to f2(At, A<fi) — f2(At). 

Fig. [5}d shows the values of variational parameters 
A(j) sa d and At sa d at the saddle points as a function of 
interaction U at half filling. One can see that the A4> sa d 
increases as U is increasing in the small U region. Af- 
ter reaching its maximum around U ~ 3.5, where both 
the difference f2(At, A<f) — il(At) and the current J be- 
come zero, A<p sa d drops sharply to zero. Note that, in 
the regions of U < 0.7 and U ~ 3.5, the difference be- 
tween Sl(At, A(j)) and f2(At) is very small, in the same 
order as the accuracy of our calculation (1.0 x 10~ 4 ). As 
in Fig. \Sjp, the absolute value of variational parameter 
At sa d is decreasing as U increases. The slope of At sa d 
curve changes abruptly at the place where the value of 
A4>sad becomes zero. 
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FIG. 3: (Color online), (a) The grand potential difference 
between Q(At, A(p) and f2(At) (denoted by □ with labeling 
on the left) and the current J corresponding to tt(At, A<f>) 
(denoted by A with labeling on the right) as a function of 
interaction U. (b) The value of variational parameters at the 
saddle points of f2(At, A(j>), A(j> sa d (* with labeling on the left 
and in unit of 7r) and At sa d ( V with labeling on the right), as 
a function of interaction U. The calculations are performed 
at half filling. 



B. Staggered flux phase in the Hubbard model 
away from half filling 

In this section, we study the stability of staggered flux 
phase when the system is away from half filling. We con- 
sider the competition between the staggered flux phase, 
superconductivity, and anti-ferromagnetism. Away from 
half-filling, the shift in the chemical potential A/i must be 
treated as variational parameter as well in order to satisfy 
the condition of thermodynamic consistency (33j j . The d- 
wave Superconducting phases [3(| H|[ can be described 
by introducing the corresponding Weiss field H' sc : 

H'sc = hsc £ ^L(c nCji + H.c), (6) 

where hsc is the strengths of the nearest-neighbor d- 
wave pairing field and Ay is the d-wave form factor. 

Let us start by comparing the superconducting and 
staggered flux states. In Fig. |4] we compare the grand 
potential of the superconducting phase ^(A/i, At, hsc) 
with that of the staggered flux phase fi(A/j,, At, A</>) for 
different values of the interaction parameter U. One can 
see that the values of the grand potential in the supercon- 
ducting state fl(A/j,, At, hsc) an d in the staggered flux 



phase f2(AyU, At, A^) are equal at half filling. This can be 
easily understood because the d-wave pairing field H' sc 
in Eq. (j6]) is connected to the staggered-flux field H' TR in 
Eq. ([31 by a particle-hole transformation at half-filling 
[43, H]. When away from half-filling, fi(A/i, At, hsc) is 
always smaller than Q (A/i, At, A<j)) for all interactions, 
implying that the superconducting state is more stable 
than the staggered flux phase in this region. For the 
large U case (U > 4.0), we can not find a VCA saddle 
point corresponding to the staggered flux phase neither 
at nor away from half filling. Here, the value of the vari- 
ational parameter A(f> is always zero as seen in Fig. |4ja) 
for U — 4.6. Fig. HJa) shows that the values of the 
grand potential tt(Afi,, At, hsc) are smaller than that of 
f2(A/x, At) and the difference between them is increasing 
as the chemical potential fi decreases. For smaller values 
of U a metastable VCA solution for staggered flux phase 
can be obtained, that is, f2(A/i, At, A(j>) is smaller than 
f2(A/i, At) with a finite value of variational parameter 
Acf>. For the intermediate U in Fig. QJb), the staggered 
flux phase disappears at some value of /i. Finally, for 
small U the metastable VCA solution for the staggered 
flux phase always exists for the whole range of chemical 
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FIG. 4: (Color online). The grand potential of super- 
conductor phase Q(A/j,, At, hsc) and staggered-flux phase 
Q(Afj,, At, A(j>), with respect to Q(Afi, At), as a function of 
chemical potential fi for various interactions U. A/i, At, Acf>, 
and hsc are variational parameters . 



potential as seen in Fig. [He) . Note that the plots in Fig. 
|2] are shown as a function of chemical potential, which 
corresponds to a doping in the range of ~ 0.1. 

In Fig. [5j the current J (a) and the value of varia- 
tional phase factor A(f> sa d (b) at the saddle points of the 
staggered flux phase are plotted as a function of dop- 
ing S for various interactions U. One can see that both 
J and Acj) sa d are zero for large U > 4.0, indicating no 
VCA solution for staggered flux phase in this region. For 
intermediate U, J and A<j) sa d decrease as the doping S 
is increasing and drop to zero at some value of doping 
(e.g., S = 0.05 for U = 3.2). For small U (e.g., U = 2.0), 
it seems like that A4> sa d keeps as a finite value (around 
0.87r) and does not drop to zero. However, the staggered 
flux phase is less stable than the superconducting phase, 
as shown in Fig. 0] 

Including the antiferromagnetic phase, we choose to 
look at the results for the values of U = 3.2, where 
the difference between Sl(Ai, A</>) and fi(At) is max- 
imal at half filling (see Fig. [3^). In Fig- El wc 
show the grand potential f2 with different variational 
parameters as a function of chemical potential //. It 
is clear that f2(A/i, At, Haf) is much smaller than both 
fi(A/z, At, h sc ) and n(A/z, At, A<j>). Thus, the antiferro- 



FIG. 5: (Color online). Current J (a) and the value of 
variational phase factor A<j> sa d (b) at the saddle points of 
fi( A/i, At, A<f>) as a function of doping S for various interac- 
tions U. 
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FIG. 6: (Color online). Grand potential with different 
variational parameters (A/i, A(j>, hsc, and /iaf) as a function 
of the chemical potential /t for U = 3.2. 



magnetic phase is the most stable state in this region of 
chemical potential. For even smaller chemical potential 
(or larger doping), the antiferromagnetic phase becomes 
unstable towards the superconducting phase. 

In summary, we find that the staggered flux phase is 
less stable than the superconducting state for all interac- 
tions and chemical potentials considered. For values of 
doping not very far away from half filling, the antiferro- 
magentic phase would be the most stable phase. 
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C. Staggered flux phase in the Hubbard model 
with nearest-neighbor interaction 

In this section, we discuss the effect of nearest- 
neighbour interaction on the staggered-flux phase in the 
Hubbard model. The nearest-neighbour interaction has 
the form: Vj2(ij) n i n ji where (ij) denotes the nearest- 
neghbour pairs. This interaction can be added to the 
Hamiltonian of Eq. ([T]). In the following, we choose 
to look at the point with a fixed on-site interaction 
U = 3.2. In Fig. [TJa), the grand potential difference 
A£l = fl(At,A(f)) — n,(At) is plotted as a function of 
nearest-neighbour interaction V at half-filling. For a pos- 
itive V, Af2 is increasing to zero as V increases, imply- 
ing that a positive repulsive V disfavors the staggered- 
flux phase. However, adding a negative V will decrease 
the value of AQ. This means that an attractive nearest- 
neighbour interaction V is favorable for the formation of 
the staggered-flux phase. However, for negative V if \V\ 
is large enough (e.g. V < —0.5), there is no stable VCA 
solution any more. The reason might be that a CDW 
order will form for very large attractive interaction and 
the staggered-flux phase is not stable any more. 



point fi = 3.6 but is still inside the Mott gap region 
(i.e, 3.5 < [i < 3.7). Afl increases rapidly to zero as 
the doping (both particle and hole) becomes nonzero 
(fi < 3.5 or [i > 3.7). For negative V = —0.5, AO 
remains nearly constant in the Mott gap region (i.e, 
—0.425 < fi < —0.375). However, interestingly, it de- 
creases as the system is dopped (i.e, fi < —0.425 or 
fi > -0.375). Note that in both Fig. Efb) and (c), the 
chemical potential corresponds to a density between 0.9 
and 1.1. Here, we do not investigate the large doping re- 
gion, in which a change of the particle sector in VCA 
is needed. Athough the grand potential difference of 
staggered-flux phase in Fig. [7jc) is decreasing with dop- 
ing, it is still hard to compete with the Antiferromagnetic 
phase. For example, the grand potential difference of 
Antiferromagnetic phase, f2(A/i, At, Haf) — f2(A/z, Ai), 
is —0.0317, which is much lower than that of staggered- 
flux phase (—0.008) at the same point. 

Therefore, our calculation shows that a negative near- 
est neighbour interaction and finite doping is energeti- 
cally favorable for the formation of staggered-flux phase. 




-0.5 -0.25 0.0 0.25 0.5 0.75 1.0 

■c- V 




FIG. 7: (Color online). (a) The grand potential differ- 
ence Q(At,A<j>) — fl(At) as a function of nearest-neighbor 
interaction V at half-filling for a fixed on-site interaction 
U = 3.2. (b) and (c) show the grand potential difference 
Q(Ajj,, At, Acfi) — Q(Afi, At) as a function of chemical poten- 
tial fi for nearest-neighbor interaction V = 0.5 and V — —0.5 
respectively. The on-site interaction is also fixed at U — 3.2. 

The effect of nearest-neighbour interaction V on sys- 
tem away from half-filling is shown in Fig. 0(b) and 
(c) for V — 0.5 and V = —0.5 respectively, where the 
Afl = £l(Afi, At, A(f>) - Q(Afi, At) is plotted as a func- 
tion of chemical potential fi. For positive V = 0.5, Af2 
increases slightly when fj, deviates from the half-filling 



D. Staggered flux phase in the three band 
Hubbard model 

Our work can be naturally extended to consider the 
three band Hubbard model to search for the existence of 
the spontaneous time-reversal symmetry breaking states. 
Such states are still under debate and have probably 
been observed in true materials 0, [45[ . In the proto- 
type realsitic systems, the high-Tc cuprates, a crucial 
role is played by the CuC>2 planes, in which the oxy- 
gen orbitals are explicit degrees of freedom. The elec- 
tron dynamics in CuC>2 plane can be described by a 
minimal model, the three band Hubbard model, contain- 
ing the copper d x 2_ y 2 orbital and oxygen p x and p y or- 
bitals. Experiments show evidence for time-reversal sym- 
metry breaking in BSCCO with photoemmision [45[ and 
in HgBa2Cu04+,5 with polarized neutron diffraction (22j. 
Early theoretical works discussed the flux phase in the 
single band model with current order and proposed it 
in conection to the pseudo-gap state in high-Tc cuprates 
[E 0, E3- I n or der to stabilize the flux phase, Varma 
shows in a mean-field approach for three band model the 
necessity to go beyond on-site interactions for an explicit 
treatment of the nearest-neighbor interactions between 
copper and oxygen 

urns. 

Orbital currents were also 
obtained by more accurate theoretical calculations for 
multi-band Hubbard models with or without the inclu- 
sion of oxygen orbitals [23|, |4(| . Within this subsection 
we show some preliminary results of VCA calculations 
for the time-reversal symmetry breaking phase in three 
band Hubbard model. 

The three band Hubbard model under consideration is 
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described by the Hamiltonian: 

H = E WW + h.c.) + J2 $(p,W + ^ c -) 

<i,j>a <j,k><7 

+^E<<+ c/ *<E n >^+ v ^ E ^ 

i j <i,j>>? 

Ha - n) E ^ + ( £ P - E ( 7 ) 

Here, dJ CT and pj^ create a hole in copper 3d x 2_ y 2 and 
oxygen 2p$ (6 — x, y) orbitals; td and e p are the on-site 
energy for Copper and Oxygen orbitals; /i is the chemical 
potential; t x L and P pp are the neighbouring Cu-0 and O- 
O hopping amplitude. Ud, U p , and Vd p are the Coulomb 
repusion when two holes sit on the same Copper orbital, 
the same Oxygen orbital, or on the neighbouring Cu and 
O orbitals. One difficulty associated with the study of 
the Cu02 planes is the large number of parameters, giv- 
ing rise to a huge phase space. Therefore, the choice of 
parameters for the three band model of copper oxides 
has been extensively discussed in the literature (47l - l49| . 
Our choice corresponds to the following realistic value 
of parameters P pp = 0.5, A = e p — = 3.0, Ud = 8.0, 
U p = 3.0, and Vd p = 0.5, where we have taken the Cu-0 
hopping t % l as unit of the energies. 
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FIG. 8: (Color online). Grand potential fi(A/i, A(f>) of the 
three band Hubbard model as a function of the values of vari- 
ational parameter A(j>. The inset shows the phase configura- 
tion used in the calculation. 

In the inset of Fig. |8] we show the Cu02 plane in 
which a 2 x 2 unit cell plaquette is taken as the reference 
system for the VCA calculations. Thus, a plaquette con- 
tains four unit cells each of them being composed of a 
Cu (d x 2_ y 2) and two O (p x , p y ) sites. The implementa- 
tion of time-reversal symmetry breaking within the three 
band model follows the description given in Sec. [TTJ In 
the present calculations we have used a reduced number 
of variational parameters, namely the shift in the chem- 
ical potential A/j, as well as the shift A<fi of the phase on 



each hopping terms, the later having the form: 
H' TR = E t%e^{dt Pja + h.c.) 

<i,j><7 

+ E 4 k p e iA ^(plp k(7 +h.c.) (8) 

<j,k>a 

A4>dp s (S — x, y) shifts the phases on the bonds between 
neighbouring Cu-0 atoms inside the 2x2 reference pla- 
quette. Similarly, A(j) pp shifts the phases on the bonds be- 
tween neighbouring 0-0 atoms whithin the same plaque- 
tte. Note that the magnitude of the shifts in the phases 
are not necessarily equal. In fact, we have performed the 
computations with four variational parameters: A.(f>d Px , 
Acj)d p , A^>pp and Afi. The direction and distribution of 
the phases on each hopping bonds are taken according 
to the "current pattern" proposed by Varma, Ref. [l7| . 
According to our calculations, none of the two current 
(phase) patterns considered in Varma's mean-field pro- 
posal [TH, [l7[ are stable. Whithin the calculations with 
four variational parameters the saddle piont of grand po- 
tential fl always locates at A(f> pp = A<fid Psc = = 0. 
Therefore no spontaneous time-reversal symmetry break- 
ing phases are found in our calculations. 

In Fig. [5] we plot the grand potential Q(A/j,, A<j>) as a 
function of Acf> for different values of /i. For the sake of 
illustration, the calculation has been performed using the 
simplifying assumption of equal phase shifts, i.e. when 
the values of A^ pp , A4>d py , and A4>d py are all equal. The 
phase (current) pattern is also shown in the inset of figure 
Fig. [5] This is the pattern denoted as the type ©// in 
the orginal work of Varma Ref. [17j . The saddle point 
in the grand potential O is obtained for zero phase on 
bonds, and turnning on phase around zero will lead to 
an increase of fl. The situation when all three A(j) are 
different is similar. 

Note that additional variational parameters such as At 
are required to complete our study on three band Hub- 
bard model, however this is out of the scope of present 
paper. 

IV. CONCLUSION 

The central aim of this paper was to extended the 
VCA method to study the time-reversal symmetry bro- 
ken phase. We proposed to use two variational parame- 
ters, the hopping At and phase factor A(j), to study the 
staggered flux phase in the 2D Hubbard model. Our cal- 
culation suggest a metastable staggered flux phase in the 
intermediate and small U regions. For the case of half 
filling, the flux phase appears to be the most stable at 
U « 3.2. This implies that the staggered flux phase, or 
at least its fluctuations, are more likely to be observed 
in the small U than in large U region. We also studied 
the stability of the staggered flux phase by comparing 
its grand potential with those of the superconductor and 
antiferromagnetic phases. We found that one of these 
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latter phases is always more stable than the staggered 
flux phase both at as well as away from half filling. 

From our results on the one-band 2D Hubabrd model 
we cannot draw a clear picture of which interactions are 
crucial to cause the orbital currents, although we clearly 
observe that a nearest-neighbour repulsion (V > 0) de- 
creases the current, that would be present in the absence 
of the intersite interaction. In contrast, an attractive 
nearest-neighbour interaction (V < 0) leads to an in- 
crease in the current comparing with the V = case. 
The staggered- flux order that emerges in the presence 
of V < was previously reported within the strong- 
coupling regime [lCf. Our numerical results demonstrate 
that such result is also applicable to the intermediate val- 
ues of U. At last, we also extented the method for time- 
reversal symmetry breaking phase to study the three 



band Hubbard model. By using only the shifts of phases 
on each hopping terms as variational parameter, we did 
not find a stable VCA saddle point for the staggered flux 
phase. 
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